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Non-Equilibrium  Phonon  Processes  and  Degradation 
in  Gigahertz  Nanoscale  Mechanical  Resonators 
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gjiafrate@ncsu.edu 

North  Carolina  State  University,  Raleigh,  NC  27695 

Report  Summary 
Objective: 

The  objective  of  this  study  is  to  identify  and  quantify  the  loss  and  degradation  mechanisms 
relevant  to  frequency  control  resonator  perfonnance  as  the  resonator  dimensions  reduce  to 
the  nanodimensional  spatial  scale.  The  Department  of  the  Army  interest  in  scaling  such 
devices  into  the  nanodimensional  region  arises  from  the  technical  notion  that  such 
resonators  can  be  designed  to  operate  at  frequencies  into  the  low  GHz  region  thus  providing 
a  low  cost,  integratable  NEMS  solid  state  device  option  for  ultrafast  frequency  control 
electronic  applications  relevant  to  jam  resistant  secure  communications.  Inherent  in  the 
nanoscale  domain  is  the  dynamics  of  non-equilibrium  phonon  processes  which  play  a 
central  role  in  orchestrating  loss  mechanisms  in  general  and  in  trade  off  of  surface  to 
volume  effects  as  the  resonator  scales  down.  In  this  study,  use  is  made  of  the  Euler- 
Bemoulli-Boltzmann  approach  to  capture  and  study  the  essence  of  the  dominant  physical 
considerations. 

Approach: 

The  technical  approach  considers  non-equilibrium  heat  generation  and  redistribution 
processes  from  mechanical  strain  during  high  frequency  NEMS  operation  beyond  the 
conventional  heat  diffusion  and  local  temperature  approximation.  A  semiclassical  phonon 
dynamical  picture  is  introduced  to  go  beyond  the  conventional  models.  Scaling  laws 
relevant  to  the  appropriate  phonon  transport  regimes  and  their  transition  boundaries  are 
delineated,  analyzed,  and  compared  with  results  of  detailed  microscopic  descriptions.  Thus, 
the  research  approach  involves  scaling  analysis,  coupled  with  analytical  and  numerical 
modeling  of  phonon  flow  and  heat  redistribution  in  nanoscale  resonators.  Specifically,  the 
following  technical  milestones  were  pursued: 

1 .  Establish  a  picture  of  the  beam  dynamics  in  a  non-equilibrium  thermodynamical 


framework.  This  requires  a  more  refined  look  at  the  approximations  of  the  Euler-Bernoulli 
equation  with  shear  corrections  to  the  beam  dynamics,  and  generalizations  of  stress-strain 
relations  in  the  presence  of  thennal  gradients. 

2.  For  intrinsic  phonon-mediated  dissipation  mechanisms,  fonnulate  scaling  laws 
based  on  relations  between  the  key  spatial  and  temporal  characteristics  of  phonons,  i.e., 

mean  free  path  1  and  phonon  relaxation  time  t ph,  and  the  relevant  physical  parameters 

of  the  NEMS  structure,  that  is  beam  thickness  t,  length  L,  and  resonator  operational 
frequency  v. 

3.  Scaling  laws  relevant  to  the  appropriate  phonon  transport  regimes  and  their 
transition  boundaries  will  be  analyzed  and  compared  with  results  of  more  microscopic 
analysis  and  theory. 

4.  Fonnulate  robust  treatment  of  phonons  under  non-equilibrium  thermodynamic 
conditions  in  flexural  nanobeams  by  introducing  Boltzmann  framework  for  description  of 
phonon  transport;  conduct  analysis  of  internal  thermal  phonon  flow  inside  the  beam. 

5.  Develop  a  reliable  quantification  of  the  Q  factor  in  presence  of  non-diflfusive 
phonon  transport. 

Accomplishments : 

•  Material  parameters  were  successfully  identified  for  the  basic  high  frequency  NEMS 
resonator  designs;  “operational  frequency  vs.  spatial  dimensions”  maps  have  been 
established  for  the  NEMS  resonator  from  the  Euler-Bernoulli  equation. 

•  Major  intrinsic  dissipative  mechanisms  related  to  thennoelastic  loss  and  phonon-phonon 
interactions  and  their  scaling  properties  have  been  identified  and  strategically  mapped 
to  provide  insightful  physical  analysis. 

•  From  the  developed  scaling  studies,  it  is  theoretically  noted  that,  in  the  1-10  GHz 
operational  frequency  regime,  the  NEMS  resonator  thennal  dynamics  routinely  goes 
beyond  the  limits  of  the  local  temperature  approximation,  not  at  all  due  to  the  often 

considered  time  constraint  of  the  high  frequency  limit  (which  requires  plp  with 

1  Ops  at  room  temperature),  but,  in  contrast,  due  to  sharp  spatial  inhomogeneity  in 
strain  pattern  induced  by  flexure  across  the  thin  (K7  with  7  «5 Onm  at  7=300  K  in 


Si)  beam  cross  section.  This  spatial  consideration  leads  to  rapid  ballistic  transfer  of 
phonons  across  the  beam  and  suppression  of  the  dissipation  mechanism  associated  with 
the  entropy  production  due  to  inter-branch  phonon-phonon  thennal  equilibration.  The 
theoretical  analysis  is  formulated  and  conducted  in  the  Boltzmann  framework  to  capture 
the  proper  phonon  dynamics. 

Journal  Publication: 

Results  of  extensive  analysis  and  conclusions  were  recently  published  by  A.  A.  Kiselev  and 
G.  J.  Iafrate  in  Phy.  Rev.  B.  77,  205436  (2008);  Publication  is  attached  as  Appendix  A; 
detail  summary,  conclusion,  and  bibliography  is  contained  therein. 
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Phonon  dynamics  and  phonon  assisted  losses  in  Euler-Bernoulli  nanobeams 

A.  A.  Kiselev* 

HRL  Laboratories  LLC,  Malibu,  California  90265,  USA 
G.  J.  Iafrate 

Department  of  Electrical  and  Computer  Engineering,  North  Carolina  State  University,  Raleigh,  North  Carolina  27695,  USA 

(Received  26  February  2008;  published  30  May  2008) 

Nonequilibrium  phonon  processes  and  related  degradation  effects  are  treated  for  a  Euler-Bernoulli  flexural 
beam  undergoing  scaling  from  a  micro  to  nanospatial  regime.  For  the  scaling  lengths  under  consideration,  the 
lowest  resonator  mode  is  in  the  frequency  range  of  1-10  GHz.  In  this  range,  it  is  found  that  the  beam  thermal 
environment  routinely  exceeds  the  limits  of  validity  for  the  local  temperature  approximation;  this  is  due  to 
sharp  inhomogeneities  in  strain  pattern  across  the  thin  beam  cross  sections  induced  by  flexural  motion  as 
opposed  to  the  often  assumed  temporal  dynamics  of  high  frequency  operation.  In  a  Euler-Bemoulli-Boltzmann 
framework,  an  analysis  of  the  internal  phonon  flow  in  the  flexural  beam  is  conducted,  and  dissipative  losses  are 
evaluated.  The  complexity  of  the  microscopic  phonon  dynamics  is  delineated  and  strategically  graphed  in 
terms  of  the  parameters  characterizing  the  flexural  beam  and  the  phonon  system  therein.  In  limiting  cases,  two 
major  intrinsic  dissipative  mechanisms  are  operative,  one  due  to  the  diffusive  spatial  redistribution  of  phonons 
resulting  in  heat  transfer  and  thermoelastic  loss,  and  the  other  due  to  thermalization  of  the  local  phonon 
population  distorted  by  strain  resulting  in  the  manifestation  of  the  Akhiezer  effect.  In  the  frequency  domain  of 
interest,  these  two  loss  mechanisms  lose  their  distinctive  character  with  decreased  spatial  scaling  and  transition 
to  a  unified  dissipative  process  governed  by  the  ballistic  phonon  transfer  across  the  beam. 

DOI:  10. 1 103/PhysRevB. 77. 205436  PACS  number(s):  66.70. -f,  63.22.Gh,  62.80.  +  f 


I.  INTRODUCTION 

Nanoelectromechanical  systems  (NEMS)  have  emerged 
as  a  very  attractive  concept  for  miniaturized  devices,  resona¬ 
tors,  and  actuators  in  application  areas  spanning  the  broad 
fields  of  metrology,  engineering,  life  science,  and 
medicine.1-2  In  particular,  microsized  resonators  and  cantile¬ 
vers  have  been  workhorse  devices  for  light-weight,  on-chip, 
and  high-quality  (high-g)  elements  for  frequency  control 
electronics  as  well  as  actuator  and  sensor  components  for 
smart  systems.  As  such  devices  are  scaled  from  the  micro  to 
the  nanolength  scale,  the  lowest  resonator  mode  of  operation 
increases  into  a  desirable  high  frequency  range  of  1-10  GHz. 
In  this  frequency  range,  the  key  question  of  concern  becomes 
whether  there  are  fundamental  intrinsic  physical  limitations 
that  come  into  play  that  might  prohibit  high -Q  operation;  this 
question  provides  the  main  motivation  for  our  work. 

Numerous  mechanisms  have  been  proposed  previously  as 
possible  limitations  to  the  quality  factors  of  micro  and 
nanoresonators.3"6  Some  mechanisms  can  be  identified  as  ex¬ 
trinsic  to  the  microresonator  and,  thus,  in  principle,  can  be 
suppressed  by  engineering  design.  For  example,  air  damping 
can  be  eliminated  in  the  vacuum-operated  resonators.  As 
well,  doping  impurities,  the  effect  of  the  attachment  and/or 
support  structure,  and  surface  absorption-desorption  pro¬ 
cesses  can  also  be  considered  as  external  mechanisms.  In¬ 
trinsic  loss  mechanisms,  on  the  other  hand,  persist  due  to  the 
inherent  material  properties  of  the  beam.  For  example,  intrin¬ 
sic  phonon-mediated  dissipation  mechanisms  are  related  to 
the  processes  of  heat  redistribution  and  entropy  production. 
This  happens  during  beam  mechanical  dynamics  because 
spatially  inhomogeneous  strain  is  induced  in  the  beam  mate¬ 
rial  by  the  motional  flexure.  Local  compression  and  dilation 
directly  affects  the  phonon  dispersion,  thereby  leading  to  the 


creation  of  nonequilibrium  phonon  distributions,  temperature 
gradients,  and,  consequently,  heat  transfer. 

In  lowest  order,  these  processes  are  most  simply  treated 
by  the  standard  theory  of  thermoelasticity  that  essentially 
relies  on  the  validity  of  the  classical  heat  diffusion  law,  the 
so-called  Fourier’s  law,5-7  8  where  heat  flow  is  proportional  to 
the  macroscopic  temperature  gradient  in  the  material  and  the 
thermal  conductivity  is  the  coefficient  of  proportionality. 
This  lowest  order  thermoelastic  methodology  does  not  in¬ 
clude  temporal  relaxation  or  higher  order  non-Fourier  heat 
flow  effects.  In  this  regard,  Cattaneo  and  Vernotte  introduced 
a  heuristic  temporal  relaxation  time  into  the  Fourier  heat 
flow  process  to  generate  a  hyperbolic  heat  equation,  which 
describes  a  nondiffusive  heat  transfer  by  means  of  the 
quickly  decaying  heat  waves.  In  particular,  with  respect  to 
the  losses  in  NEMS,  Sun  et  al.9  built  on  the  temporal  relax¬ 
ation  approach  by  incorporating  a  non-Fourier  heat  flow  term 
into  their  analysis;  it  is  worth  mentioning  that  for  the  mag¬ 
nitude  of  parameter  values  numerically  considered  in  their 
paper,  the  effect  of  non-Fourier  heat  transfer  is  minimal.  Fur¬ 
thermore,  in  the  work  of  Joshi  and  Majumdar,10  the  “Fouri¬ 
er’s  law”  and  “Cattaneo”  nonstationary  solutions  for  a  simple 
one -dimensional  geometry  were  directly  compared  to  the 
analysis  from  a  microscopic  Boltzmann  treatment;  it  was 
comparatively  noted  that  both  the  Fourier  and  Cattaneo  re¬ 
sults  were  inadequate  on  short  spatial  and  temporal  scales. 
This  comparative  analysis  points  to  the  notion  that  micro¬ 
scopic  phonon  dynamics  should  be  accounted  for  holistically 
in  order  to  facilitate  a  more  complete  and  accurate  analysis 
of  the  heat  redistribution  processes  in  nanodimensional  struc¬ 
tures.  It  is  worth  mentioning  that  the  overwhelming  com¬ 
plexity  of  solving  directly  the  Boltzmann  equation,  even  for 
relatively  simple  configurations,  brought  to  existence  the  ad¬ 
vanced  simplifying  schemes  to  numerically  treat  heat  redis- 
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tribution  processes,  such  as  the  nonlocal  heat  conduction 
model,11  and  the  ballistic-diffusive  heat-conduction 
equations.12 

Another  intrinsic  phonon-mediated  loss  mechanism  is 
known  as  the  Akhiezer  effect.13  This  effect  is  most  com¬ 
monly  considered  in  relation  to  the  attenuation  losses  suf¬ 
fered  by  ultrasonic  waves  in  dielectric  crystals  at  elevated 
temperatures.  In  essence,  due  to  anharmonicity,  a  strain 
modulates  the  phonon  frequencies  (and,  most  important,  dif¬ 
ferently  for  different  phonon  modes);  thus,  the  original  pho¬ 
non  distribution  becomes  distorted  and  requires  a  micro¬ 
scopic  time  to  reestablish  phonon  equilibrium  locally.  This 
irreversible  process  results  in  an  absorption  of  elastic  energy 
to  generate  entropy  in  the  phonon  subsystem.  A  simple  quali¬ 
tative  treatment  of  this  process  was  reported  by  Bommel  and 
Dransfeld.14  A  more  advanced  microscopic  treatment  was  de¬ 
veloped  by  Woodruff  and  Ehrenreich.15 

Unlike  ultrasonic  waves  in  infinite  solids,  where  the  scale 
of  the  spatial  inhomogeneity  of  the  strain  pattern  is  defined 
solely  by  the  plane-wave  frequency  through  the  bulk  disper¬ 
sion  law,  the  strain  pattern  associated  with  flexural  waves  in 
thin  beams  and  membranes  is  aptly  influenced  by  the  bound¬ 
ary  conditions  associated  with  the  mechanical  motion  of  the 
finite  material  domain.  As  a  result,  the  spatial  and  temporal 
scales  of  the  strain  pattern  can  be  controlled  more  indepen¬ 
dently.  Specifically,  sharp  spatial  inhomogeneity  in  strain 
pattern  on  the  relevant  scale  of  the  thermal  phonon  mean  free 
path  can  be  routinely  induced  by  the  flexure  in  the  cross 
section  of  the  gigahertz  thin  beam.  This  simple  observation 
has  direct  and  important  consequences  for  the  magnitude  and 
variation  of  intrinsic  losses  in  nanoscale  resonators. 

In  the  general  framework  of  this  paper,  the  detailed  pho¬ 
non  dynamics  is  established  for  Euler-Bernoulli  flexural 
beams;  as  well,  phonon-related  losses  and  degradation 
mechanisms  are  identified  for  a  mechanical  resonator  beam 
configuration  as  the  beam  size  reduces  to  the  nanodimen¬ 
sional  spatial  scale  and  the  frequency  of  operation  rises  into 
the  low  GHz  spectral  region.  The  technical  approach  consid¬ 
ers  nonequilibrium  heat  generation  and  redistribution  pro¬ 
cesses  from  mechanical  strain  during  high  frequency  NEMS 
operation  beyond  the  conventional  heat  diffusion  and  local 
temperature  approximation.  Specifically,  a  semiclassical 
phonon  dynamical  picture  is  introduced,  allowing  detailed 
microscopic  analytical  and  numerical  modeling  of  phonon 
flow  and  heat  redistribution  in  Euler-Bernoulli  nanoscale 
resonators.  Relevant  mechanistic  phonon  transport  regimes 
and  related  transition  boundaries  are  identified  by  delineating 
the  appropriate  scaling  laws. 

In  Sec.  II,  the  salient  features  relevant  to  the  dynamics  of 
the  Euler-Bernoulli  flexural  beams  are  reviewed  and  specific 
attention  is  devoted  to  the  structure  of  the  spatially  inhomo¬ 
geneous  time  dependent  strain  pattern  developed  inside  beam 
material.  Section  III  is  focused  on  the  formulation  and  analy¬ 
sis  of  phonon  dynamics  under  the  nonequilibrium  thermody¬ 
namic  conditions  in  flexural  nanobeams  within  a  Boltzmann 
framework;  in  particular,  a  reliable  quantification  of  the  nan¬ 
oresonator  Q  factor  is  developed  in  the  presence  of  nondif- 
fusive  phonon  transport.  In  Sec.  IV,  a  numerical  analysis  of 
the  phonon  assisted  losses  in  the  double  clamped  flexural 
NEMS  resonator  is  provided  as  an  illustrative  example  of  the 


methodology.  The  functional  dependence  of  the  loss  on  the 
characteristics  of  the  flexural  and  phonon  systems  is  estab¬ 
lished  and  graphically  portrayed.  Specifically,  domains  of 
phonon-related  dissipative  mechanisms713  are  identified  and 
analyzed  in  terms  of  their  scaling  attributes  relative  to  the 
key  spatial  and  temporal  characteristics  of  the  phonons,  i.e., 
mean  free  path  and  phonon  relaxation  time,  as  well  as  the 
relevant  physical  parameters  of  the  NEMS  structure,  i.e., 
beam  thickness,  length,  and  resonator  operational  frequency. 
In  Sec.  V,  results  are  summarized,  and  further  considerations 
concerning  the  integrity  of  Q  factor  are  discussed. 


II.  PICTURE  OF  THE  BEAM  DYNAMICS 


Flexural  modes  available  to  the  double  clamped  resonator, 
cantilever  (for  large  length  to  thickness  aspect  ratios  L/t 
>10),  or  in  the  model  of  an  infinite  beam  (then  L  is  the 
period  of  the  mode)  can  be  described  analytically  by  the 
Euler-Bernoulli  equation  of  solid  mechanics  in  the  limit  of 
no  energy  dissipation  and  in  the  continuum 
approximation.1,16  Applying  this  long-beam  approximation 
of  the  Euler-Bernoulli  theory,  one  arrives  at  the  equation  for 
the  transverse  y  displacement  u  (z ,  t)  for  an  undriven  beam 
along  the  direction  of  its  length  z  as 


cfu  d1  cfu 

pA  t  -j  0^1  2  ~  0- 

dr  dz  dz 


Here,  p  is  the  material  density,  A  is  the  area  of  the  beam 
cross  section,  E  is  the  Young  modulus,  I  is  the  cross- 
sectional  area  of  the  moment  of  inertia  (bending  moment  of 
inertia),  and  El  is  known  as  the  flexural  rigidity  of  the  beam. 
The  boundary  conditions  specify  the  beam  clamping  condi¬ 
tions  for  finite  structures.  In  particular,  for  the  rigidly 
clamped  boundary  condition  u= 0  and  u'  =  duldz= 0.  Thus, 
for  the  double  clamped  resonator  (clamped  at  beam  ends  at 
z  =  0  and  z=L),  the  eigensolution  for  a  homogeneous  beam  of 
a  rectangular  cross  section  wt  (i.e.,  A  =  wt  and  /=wf3/12), 
can  be  written  in  phasor  form  as 

un(z,t)  =  Un(z)ein'f 

where 


Un(z)  =  C^cos^z)  -  cosh(£„z)]  +  C)!2[sin(fc„z)  -  sinh(£„z)], 
with  frequency 


Here,  kn  is  chosen  to  satisfy  the  boundary  conditions  and  is  a 
solution  of  the  transcendental  equation, 

cos(knL)cosh(knL)  =  1. 

It  is  found  that  knL~  4.730,  7.853,  and  10.996  for  the  first 
three  modes.  The  fundamental  frequency  11,  of  the  oscillator 
rises  as  the  resonator  size  and/or  aspect  ratio  (length  to  thick¬ 
ness)  is  reduced  as  desired. 

It  is  noted  that  while  the  Euler-Bernoulli  theory  works 
well  for  the  relatively  long  beams,  the  importance  of  the 
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Timoshenko  rotational  inertia  and  shear  corrections17  to  the 
beam  dynamics  increases  progressively  for  smaller  Lit  ra¬ 
tios.  Also,  our  treatment  of  double  clamped  beams  here  is 
limited  to  small  flexural  vibrations,  which  assumes  that 
|t/(z)|<f;  nonlinearities  are  set  into  play  through  the  qua¬ 
dratic  dependence  on  U  in  the  beam  lengthening  term,1  A L 
~  f idz[U' (z)]2/2,  as  larger  flexural  vibrations  affect  the 
resonator  dynamics.  For  particular  classes  of  flexural  sys¬ 
tems,  such  as  carbon  nanotubes  with  naturally  small  thick¬ 
nesses,  even  very  moderate  displacements  would  lead  to 
resonator  dynamics  beyond  the  linear  regime.18  On  the  other 
hand,  for  cantilevers,  the  significance  of  this  deviation  from 
linearity  is  largely  alleviated. 

The  structure  of  the  strain  field  corresponding  to  a  par¬ 
ticular  flexural  mode  can  be  established  as  follows:  For  thin 
beams  in  the  linear  regime,  the  situation  is  very  simple,16  the 
only  nonzero  components  of  the  strain  tensor  are  defined  by 
the  mode  deflection  amplitude  u(z)  as 

o*2 u 

Uzz  =  ~y^2,  (2) 

oz 

and 

It XX  ~  ttyy  —  ~  tZU^z,  (3) 

where  cr  is  known  as  the  Poisson  ratio.  The  trace  of  the  strain 
tensor  is  then, 

2  «,,  =  —  (  1  -  2a)yu" .  (4) 

i 

The  linear  energy  density  along  the  beam  is  given  in  pha- 
sor  form  by 

hJz)  =  n^[U(z)f+f[U"(z)]\ 

where  the  first  term  is  the  kinetic  energy  and  the  second  term 
is  the  potential  energy  of  flexural  deformation  of  the  beam. 
The  kinetic  energy  density,  in  phasor  form,  is  proportional  to 
the  square  of  the  fundamental  mode  eigenfunction  l]\  and 
the  flexural  deformation  energy  is  proportional  to  [I/"].2 
Both  components  of  the  energy  distribution  are  highly 
peaked  along  the  beam  and  scale  in  accordance  with  L.  How¬ 
ever,  as  L  is  reduced,  the  relative  amounts  of  the  kinetic  and 
flexural  energies  along  the  beam  remain  unchanged  for  the 
Euler-Bernoulli  normalized  eigenmode  Ui ;  thus,  the  kinetic- 
flexural  energy  partition  could  vary  with  L  only  by  consid¬ 
ering  the  higher  order  corrections  to  the  Euler-Bernoulli 
equation. 

Finally,  the  total  energy  of  the  mechanical  oscillation  is 
given  by 

W=  |  dzhw(z).  (5) 

J  L 

III.  PHONONS  UNDER  NONEQUILIBRIUM 
THERMODYNAMIC  CONDITIONS 
IN  FLEXURAL  NANOBEAMS 

A  quantitative  description  for  the  processes  of  heat  gen¬ 
eration  and  redistribution  in  thin  flexural  beams  is  developed 
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by  applying  the  Boltzmann  formalism  to  capture  the  dynam¬ 
ics  of  the  nonequilibrium  phonon  distribution  resulting  from 
the  spatially  inhomogeneous  mechanical  compression  and 
dilation.  The  appropriate  multidimensional  integro- 
differential  equation  is  derived  for  the  phonon  distribution 
function,  including  the  modulation  of  phonon  frequencies  by 
spatially  inhomogeneous  strain,  phonon  ballistic  transfer,  and 
scattering  processes,  leading  to  the  thermalization  of  the  pho¬ 
non  ensemble.  This  treatment  allows  for  a  reliable  quantifi¬ 
cation  of  the  phonon  flow  in  the  resonator,  and  an  analysis  of 
the  quality  factor  Q  in  the  presence  of  nondiffusive  phonon 
transport. 


A.  Phonon  Hamiltonian  in  the  presence  of  strain 

The  space  and  time-modulated  Hamiltonian  for  a  phonon 
with  the  wave  vector  q  can  be  written  as 

H=  fico(q,r,t)  =  fuo0(q)[l  +  a(q,r,f)]; 

here,  for  the  simplicity  of  the  notation,  q  includes  both  the 
wave  vector  and  index  of  the  phonon  branch.  In  the  presence 
of  strain,  the  lowest  order  expression  for  the  modulation  co¬ 
efficient  a  is  given13  by 


a(q,r,f)  =  -  2  y/*(q)«,*(iM),  (6) 

ik 

where  ujk  is  the  strain  tensor  and  yiJt(q)  is  the  generalized 
Griineisen  tensor ,15  The  only  nonzero  components  of  the 
strain  tensor  in  the  beam  are  those  in  Eqs.  (2)  and  (3)  above. 
Then,  one  can  write  explicitly  from  Eq.  (6), 


a(q,y,z,t)  ■■ 


(i  +  <j)y7,(q)  ~  y„(q) 


yu"(z,t).  (7) 


In  replacing  all  of  the  diagonal  tensor  components  y„  with  a 
single  Griineisen  constant  y,  we  can  reduce  Eq.  (7)  further, 
so  that, 


a  =  -  y2  Uu  =  (1  -  2 o)  yyu"{z ,  t ) .  (8) 


The  simplification  that  y  assumes  a  scalar  and  ^-independent 
form  was  invoked  historically  in  the  early  development  of 
the  bulk  thermal  expansion  theory  and  has  served  well  in 
simple  situations  requiring  only  the  averages  over  the  bulk 
phonon  spectrum;  however,  in  the  limiting  case  of  a  nanore¬ 
gime,  the  averaging  with  respect  to  A,y„  and  yzz  is  not  easily 
justified,  especially  when  one  needs  to  accommodate  the  in¬ 
herent  physically  present  spatial  anisotropy  of  the  system  in 
the  analysis  of  the  degradation  processes.  Moreover,  it  fol¬ 
lows  that  in  choosing  yiy( q)  as  a  single  constant,  it  artificially 
precludes  the  manifestation  of  the  Akhiezer  effect;  in  order 
to  preserve  the  physical  presence  of  the  Akhiezer  phenom¬ 
ena,  potentially  important  at  the  NEMS  spatial  and  frequency 
scale  yet  keep  the  analysis  tractable,  one  could  allow  for  two 
different  constants  y^  and  y2  to  mimic  the  variation  in  tensor 
yl;(q) — an  approach  originally  proposed  by  Bommel  and 
Dransfeld.14  This  two  parameter  approach  is  utilized  later  in 
this  section. 
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B.  Boltzmann  equation 

The  quantitative  description  of  the  heat  generation  and 
redistribution  processes  in  thin  flexural  beams  is  developed 
using  the  Boltzmann  formalism.  The  linearized  Boltzmann 
equation  is  derived  to  capture  the  responsive  development  of 
the  nonequilibrium  phonon  distribution  resulting  from  the 
spatially  inhomogeneous  mechanical  compression  and  dila¬ 
tion  during  flexure.  For  the  instantaneous  distribution  func¬ 
tion  N(q,r  ,t)  for  the  population  of  the  phonons  in  the  beam, 
the  Boltzmann  equation  is  expressed  as 

/  dN\  dN  dN  d(0  dN  da> 

\  dt  /coll  dt  dr  dq  dq  dr 


Splitting  N  into  the  “equilibrium”  part  N0(fuol  kT0),  corre¬ 
sponding  to  the  local  phonon  frequencies  co  defined  by  the 
spatially  inhomogeneous  strain  and  the  true  equilibrium  tem¬ 
perature  T0,  and  the  nonequilibrium  part, 

N1(q,r,t)=N(q,r,t)-N0(h(o/kT0),  (9) 


one  can  linearize  the  original  Boltzmann  equation  for  the 
phonon  distribution  function  to  arrive  at19 


/  dN\  dNi  ,hu>r\da  dN\ 

—  = - l-+N' — +  v0- - L. 

\df/coll  dt  0  kT0  dt  dr 


(10) 


Here, 

Nq(u>,T)  =  No{hu>lkT)  =  [exp (hw/kT)  -  1]_1 


is  the  thermal  equilibrium  distribution  function,  N'Q  is  its  de¬ 
rivative  in  respect  to  x=hu>/kT,  and  v0=d&>0/dq  is  the  mode 
velocity. 

Invoking  the  effective  relaxation  time  approximation,  one 
can  write  for  the  collisional  term. 


along  the  y  axis  in  the  plane  of  the  beam  cross  section.  Lin¬ 
earizing  the  collision  term  of  Eq.  (11)  with  help  of  Eq.  (9), 
and  then  seeking  oscillatory  solutions  of  Eq.  (10)  for  N\  with 
flexural  fundamental  frequency  (1,  we  obtain  for  N\, 


(l-iftrph)Ali-u07-ph  cos 


6N\  =  K— 

1  0  kTn 


iflTpha- 


(12) 


Here,  the  spatial  variation  of  a  is  given  by  Eq.  (8)  and  the 
temperature  difference  7)  =  7’-  70  is  defined  by  the  total  en¬ 
ergy  balance  in  the  scattering  processes,  i.e., 

„  /  dN\ 

dqhaA —  =0.  13) 

\  dt  /coll 


Further  analytical  progress  can  be  made  by  making  several 
simplifying  assumptions:  (i)  assume  the  Debye  model 
Cfj{](q)  =  c()q  with  the  ^-independent  velocity  u0  =  c0  for  the 
phonon  spectrum;  and  (ii)  assume  that  the  relaxation  time  rph 
is  also  phonon  mode-independent.  Then,  in  treating  the  col¬ 
lisional  term  in  the  relaxation  approximation  and  expanding 
Eq.  (13)  for  the  energy  balance  to  first  order,  one  obtains. 


1 


J  d  cos  6  J 


dqq2fia>{q)N  \ 


2 

hm0(q) 

1  dqcf 

.  kT0  _ 

kNn. 


(14) 


The  equation  for  the  right  hand  side  (rhs)  is  easily  recog¬ 
nized;  the  specific  heat  of  phonon  mode  q  is 


ldN\  N-Nq(cj,T) 

(ii) 

S(q)  =  - 

hco0(q) 

\  dt  J  coj]  Tph 

[  kT0  \ 

No¬ 


where  rph(q)  is  the  phonon  relaxation  time.  There  are  a  num¬ 
ber  of  subtle  issues  associated  with  Eq.  (11).  First,  it  is  a 
heuristic  approximation  to  the  realistic  microscopic  descrip¬ 
tion  for  the  phonon-phonon  interactions.20  Even  then,  actu¬ 
ally  two  classes  of  phonon-phonon  scattering  processes  take 
place:  The  normal  processes  that  conserve  total  momentum 
of  the  colliding  pair  and  umklapp  processes  that  do  not  con¬ 
serve  total  momentum.  As  written  in  Eq.  (11),  it  is  implicitly 
assumed  that  the  umklapp  processes  dominate  so  that  the 
asymptotic  phonon  distribution  can  be  characterized  by  a 
single  parameter — local  equilibrium  temperature  T.  Further¬ 
more,  this  local  equilibrium  temperature  T  should  not  be  de¬ 
fined  by  averaging  the  phonon  distribution  only  at  a  single 
spatial  point;  instead,  it  should  be  obtained  by  averaging 
over  the  area  of  the  size  of  about  /ph.  Thus,  it  is  clear  that  T 
should  approach  TQ  as  the  beam  thickness  becomes  smaller 
than  /ph.  Nevertheless,  it  was  shown15  that  treating  T  as  a 
single-point  quantity  should  produce  comparable  results  for 
the  losses  even  in  the  limit  of  r</ph. 

Now,  we  consider  a  case  of  substantially  different  spatial 
scales  L>t,  suitable  for  the  Euler-Bemoulli-Boltzmann  reso¬ 
nator.  In  particular,  this  allows  for  the  treatment  of  all  ther¬ 
mal  transfer  as  essentially  one-dimensional — taking  place 


and  the  heat  capacity  of  the  phonon  subsystem  (per  unit  vol¬ 
ume)  is 


C=^J  dqq2S(q); 

therefore,  the  rhs  of  Eq.  (14)  reduces  to  CTl. 

Furthermore,  assuming  the  Griineisen  tensor  y  to  be  a 
^-independent  scalar  constant  as  given  by  Eq.  (8),  one  can 
multiply  Eq.  (12)  by  a  phonon  energy  and  then  integrate  over 
modulus  q  to  get  an  equation  for  the  partially  integrated  pho¬ 
non  distribution  function. 


f{y,  cos  9,z) 


dqq2ha>oN\. 


Explicitly,  one  gets, 

(1  -  iilTph)f-  lph  cos  0f'y  =  C[T |  -  i(lTphT0a(y,z)\, 

(15) 


with  Eq.  (14)  reducing  to. 
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Using  Eq.  (8)  to  express  a(y,z )  via  y  and  the  flexural  mode 
deflection  U{z),  it  is  convenient  to  further  introduce  a  func¬ 
tion  f(y,  cos  0)  as 

/O', cos  8,z)  =  cr0(l  -2o)tU"{z)f{y, cos  6), 

where  y=ylt.  Here,  the  advantage  lies  in  the  explicit  factor¬ 
ization  of  the  z  dependence,  i.e.,  /  is  functionally  the  same 
along  the  beam  length  and  satisfies  the  equation. 


(1  -  iflTph)/-  7ph  COS  8fy  =  F-  iO,Tphyy,  (16) 

where  Zph=  /ph/ ?  is  a  dimensionless  system  parameter  and 
2F(y)  =  Jd  cos  Of.  Note  that,  alternatively,  one  could  decide 
to  introduce  y=y/Zph,  T=t/L h,  and  /ph=  1,  but  our  choice 
seems  to  result  in  a  slightly  more  transparent  scaling  law 
analysis  for  loss  in  the  resonator.  Lastly,  Eq.  (16)  needs  to  be 
accompanied  by  physically  appropriate  boundary 
conditions21  and  then  solved  either  numerically  directly  or  as 
a  series  expansion. 

In  addressing  the  boundary  condition,  it  is  instructive  to 
present  /  as  a  sum  of  two  functions:  one  symmetric  with 
respect  to  the  angular  argument  of  cos  6  and  the  other  anti¬ 
symmetric  with  respect  to  cos  O',  therefore. 


f  =  fs+fa- 

Then  2 F=  Jd  cos  0f=  Jd  cos  0fs  alone,  whereas  fa  can  be 
expressed  via  fs  with  an  equation  for  fs  given  by, 

(1  -  in,Tpb)2fs-?ph  cos2  0f'  =  (\  -  iilTph)[F -  aiTpbyy]. 

Thus,  in  the  case  of  perfectly  reflecting  boundary 
conditions,21  function  fa  vanishes  at  the  beam  edges,  which 
corresponds  to/s'=0  at  y  =  ±  1/2.  Expressing  the  solution  for 
fs  as  an  expansion, 

fs(y,  0)  =  fs,m(0)sin(2m  +  1)tt y,  (17) 

m 

and  noticing  that,  identically, 

y  =  2  ym  sin(2 m  +  l)7ry, 

m 

with  coefficients  ym  =  2(-l)m/  Ti2{2m  + 1)2,  it  follows  that  the 
coefficients  fs  m  and  their  related  angular  integrals  Fm  can 
then  be  calculated  analytically  without  much  difficulty. 


C.  Phonon  related  losses  in  flexural  beams — a  single 
constant  scalar  Griineisen  coefficient 

As  already  mentioned  in  previous  discussion,  the  problem 
of  losses  in  flexural  beams  is  conceptually  similar  to  a  prob¬ 
lem  of  sound  attenuation  in  bulk  dielectric  crystals.14,15  In 
the  bulk  discussions,  an  infinite  domain  with  solutions  in  the 
form  of  plane  waves  was  analyzed;  in  this  paper,  the  finite 
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boundaries  force  solutions  on  the  finite  domain  in  the  pres¬ 
ence  of  spatial  inhomogeneities.  The  energy  conservation 
law  enforces  the  condition  that  the  average  rate  at  which  the 
energy  is  removed  from  the  flexural  vibration  is  equal,  in  the 
steady  state  process,  to  the  rate  of  transferring  energy  from 
the  phonon  system  to  the  heat  bath22  through 


z=-2 


(18) 


Using  Eq.  (10)  to  express  the  collisional  term  in  Eq.  (18)  via 
the  time  derivative  of  N  and  drift  and  diffusion  terms,  and 
further  removing  terms  that  time  average  to  zero  in  the 
steady  process,  it  is  found  that. 


With  inverse  resonator  quality  factor  Q  formally  defined  as 
the  fraction  of  the  oscillation  energy  lost  per  radian  of  flex¬ 
ural  vibration,  i.e.. 


Q~'  = 


z 

aw’ 


(20) 


the  final  result  for  quality  factor  is  found  after  substituting 
the  expression  for  the  total  flexural  energy  W  given  by  Eq. 
(5)  and  the  average  rate  of  energy  transfer  Z  defined  by  Eq. 
(19)  into  Eq.  (20). 

In  utilizing  the  previously  made  approximations:  (i)  a 
large  beam  aspect  ratio,  (ii)  the  Debye  model,  and  (iii)  a 
single  constant  scalar  Griineisen  coefficient,  one  can  perform 
a  number  of  partial  integrations  to  reduce  Eq.  (19)  to 

Z=a(1~2<7)2Cr(,yVl’f  J  dz[U"(z)]2J  dyyIm{F(y)j, 

(21) 

in  terms  of  function  F  for  a  beam  mode  U.  This  gives 

QTl  =  6(1  -  2cr)2CT0E~1yJ  dyy  lm{F(y)j.  (22) 

With  Ni  (or  /  and  its  cos  0-integral  F)  determined,  phonon- 
related  losses  in  the  beam  can  be  numerically  evaluated.  In 
particular,  if  /  is  expressed  as  a  series  expansion  of  the  form 
given  in  Eq.  (17),  then, 

QTl  =  6(1  -  2cr)2CT0E~l y^j  ym  Irn{Em}.  (23) 


Here,  it  is  stressed  that  when  the  assumption  of  one¬ 
dimensional  thermal  transport  (case  of  a  relatively  long 
beam,  that  is  L>t)  is  not  valid,  Nl  cannot  be  factorized  to 
effectively  separate  the  coordinate  z  along  the  beam,  thus, 
leaving  us  with  a  more  complex  system  of  partial  differential 
equations  to  be  solved  numerically. 

D.  Losses  in  the  presence  of  the  Akhiezer  effect 

Adapting  a  heuristic  based  on  only  a  single  Griineisen 
constant  precludes  the  manifestation  of  the  Akhiezer  effect. 
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Therefore,  a  heuristic  which  allows  us  to  capture  the 
Akhiezer  effect  is  now  developed.  As  a  starting  point  in  con¬ 
sidering  losses  in  acoustic  waves  in  solids,  Bommel  and 
Darnsfeld14  have  split  all  phonon  modes  into  two  distin¬ 
guished  groups — those  experiencing  positive  temperature 
changes  and  those  experiencing  negative  or  no  temperature 
changes.  Following  Ref.  14,  we  too  assumed  that  all  phonon 
states  can  be  split  into  two  global  groups  labeled  “1”  and 
“2,”  with  portion  pl  of  all  phonons  belonging  to  the  first 
group  and  portion  p2=l~P\  belonging  to  the  second  one. 
Both  groups  are  further  assumed  to  have  an  exactly  identical 
spectrum  in  absence  of  strain.  For  simplicity,  the  effect  of  the 
applied  strain  will  be  controlled  by  mode-independent  (for 
each  of  two  groups)  Griineisen  constants  j\  and  y2.  Then 
(the  symmetric  portions  of)  the  phonon  distribution  functions 
fis  in  each  group  (/=  1,2)  can  be  found  as  a  solution  to  the 
system  of  equations, 

(1  -iiiTph)2fls-?ph  cos2  0f;s  =  {\-iD,Tph)[F-iilTphy^], 
with  cos  0-integrals  defined  as 

2Fl.j„cosV„,  F  =  plF!+p2F2. 

Solving  this  system  analytically  via  series  expansion,  allows 
to  obtain  Q  as 

Q-1=6(l-2a)2CT0E-1^ym  lm{ply1Flm+ p2y2F2m}. 

m 

(24) 


IV.  NUMERICAL  ANALYSIS  AND  DISCUSSION 

In  Fig.  1,  graphical  analysis  for  the  phonon-related  losses 
(2-1)  is  presented  in  frequency-beam  thickness  (I2,f)  coor¬ 
dinates  (logarithmic  scales)  for  the  double  clamped  flexural 
resonator.  The  magnitude  of  the  loss  is  indicated  by  the  color 
code:  light  for  high  Q  or  low  losses;  dark  for  low  Q  or  high 
losses.  The  thickness  and  the  frequency  axes  to  the  left  and 
bottom  of  the  graph  are  scaled  in  dimensionless  units,  f//ph 
and  Ilrph,  respectively.  Corresponding  values  for  the  right 
and  top  axes  are  scaled  specifically  for  crystalline  silicon  (Si) 
at  room  temperature.  Taking  for  Si  the  standard  values  of  the 
sound  velocity  c  =  5.86  km/s,  the  thermal  conductivity  k 
=  160  W/mK,  and  the  heat  capacity  (per  unit  volume)  C 
=  1.64X  106  J/m3  K,  one  obtains  for  the  characteristic  pho¬ 
non  microscopic  time  constant  rph  —  1 0  ps,  which  corre¬ 
sponds  to  the  phonon  mean  free  path  Zph=cTph  — 50  nm. 
Other  Si-specific  numbers  used  in  the  present  analysis  are  as 
follows:6  density  p= 2330  kg/m3,  Young’s  modulus  E 
=  1.69X10U  kg/m  s2,  and  the  (linear)  thermal  expansion 
coefficient  aT=  2.8  X  10-6  K_l.  These  define  the  value  of  the 
average  Griineisen  constant  y  as  —0.29.  Also,  the  Poisson 
ratio  cr=0.27  is  used. 

In  these  graphs,  the  position  of  beams  with  length-to- 
thickness  aspect  ratios  of  LI  t  =  5,  and  10  is  marked  by  dotted 
lines  as  indicated.  As  appropriate  to  the  Euler-Bernoulli  ap¬ 
proximation,  the  analysis  of  the  thermal  phonon  flow  is  con¬ 


ducted  for  the  case  L  9>  t  as  well  as  L  9>  Zph;  in  this  case,  the 
simplifying  assumption  is  invoked  that  phonon  heat  transfer 
is  primarily  one-dimensional.  Thus,  only  beam  sizes  with  L 
S5 1  (and  with  t  not  too  much  less  than  f=/ph)  are  treated 
rigorously  (lower  left  part  of  the  graphs).  Serious  deviation 
can  be  expected  when  these  assumptions  are  not  satisfied. 

The  Q  values  in  Fig.  1  are  calculated  from  Eq.  (24),  which 
is  generally  parametrized  to  reflect  losses  in  the  presence  of 
the  Akhiezer  effect.  Specifically,  three  distinguishable  cases 
are  considered  in  the  framework  of  the  generalized  Bommel- 
Dransfeld  approach:  (a)  y\  =  y2=  y,  (b)  y !=  y=-y2,  and  (c) 
7i  =  2yA 0,  y2=0.  For  all  cases,  we  have  chosen  the  propor¬ 
tions  of  the  two  groups  of  phonons  to  be  P\=p2  =  0.5. 

On  physical  grounds,  results  reflected  in  these  graphs  can 
be  conveniently  analyzed  with  the  help  of  the  well-known 
Zener  model.3  As  a  background  discussion  relevant  to  this 
model,  macroscopic  relations,  such  as  the  linear  stress-strain 
relation,  with  Young  modulus  E  serving  as  a  coefficient  of 
proportionality  often  deemphasize,  for  simplicity,  the  aspect 
of  temporal  evolution.  More  generally,  £  is  a  function  of  the 
modulation  frequency  12  and  is  neither  necessarily  constant, 
nor  necessarily  real,  thus,  giving  rise  to  loss  and  memory 
effects.  The  reason  for  this  temporal  aspect  is  that  there  is  a 
finite  time  necessary  for  the  statistical  correlations  in  the  bath 
to  adjust  to  the  system  motion;  therefore,  friction  (internal 
friction)  arises  as  a  result  of  the  nonlinear  interaction  of  the 
system  with  the  bath,  leading  to  transfer  of  excess  potential 
and  kinetic  energy  from  system  to  the  bath.  The  magnitude 
of  friction  is  typically  proportional  to  the  ensemble  average 
of  the  system  velocity  (e.g.,  duldt). 

A  simple  (and  popular)  phenomenological  generalization 
of  the  static  Hooke’s  law  allowing  build  up  time  for  strain 
and  stress4  leads  to  Lorentzian  peaks  in  both  the  imaginary 
part  of  Young’s  modulus  £(I2)  and  loss  factor  Q~x .  This  gen¬ 
eralization,  which  is  the  basis  of  Zener’s  model,  suggests  that 
in  the  presence  of  the  relaxation  processes  with  characteristic 
time  t,  the  inverse  quality  factor  (losses)  as  a  function  of  the 
oscillation  frequency  12  can  be  described  by  a  universal  de¬ 
pendence. 


It  is  clear  that  this  functional  dependence  has  a  peak  at 
12peak=  t-1  .  On  a  log-log  scale,  this  peak  manifests  itself  as  a 
distinctive  symmetric  hump  with  45°  slopes:  for  12 r<S  1,  the 
processes  are  “isothermal”  with  a  high  level  of  environmen¬ 
tal  adjustment  achieved  during  operation  while  in  the  oppo¬ 
site  limit  of  12t5>1,  the  motion  is  adiabatic  as  any  adjust¬ 
ment  in  the  bath  is  too  slow  to  result  in  notable  friction.  In 
the  presence  of  several  relaxation  mechanisms,  several  peaks 
of  different  magnitudes  (defined  by  the  value  of  A)  will  be 
observed. 

For  the  fixed  beam  thickness,  case  (a)  demonstrates  a  well 
defined  single  peak  in  the  losses  as  a  function  of  frequency. 
On  the  log-log  (12 ,  t )  map,  peak  position  follows  a  piecewise 
linear  dependence  with  a  crossover  region  corresponding  ap¬ 
proximately  to  r=/ph.  At  large  beam  thicknesses,  the  slope  of 
the  linear  dependence  corresponds  to  fpeak  'x  I2pl!.^  while  at 
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FIG.  1.  Magnitude  of  phonon  assisted  losses  in  Si  double 
clamped  resonator  at  room  temperature  as  a  function  of  the  flexural 
frequency  fl  and  the  beam  thickness  t.  Color  code:  light — high  Q , 
low  losses;  dark — low  Q,  high  losses,  (a)  single  Griineisen 
constant — pure  case  of  thermoelastic  losses.  At  V’  the  character 
of  phonon  transport  changes  from  diffusive  to  ballistic  (Fourier  vs 
non-Fourier  heat  transfer),  (b)  Two  exactly  equivalent  groups  of 
phonon  states  with  opposite  Griineisen  constants — losses  dominated 
by  the  Akhiezer  effect,  (c)  Interplay  of  both  mechanisms  (here,  we 
assume  equality  of  the  thermoelastic  and  microscopic  phonon- 
phonon  prefactors).  Strictly  speaking,  only  the  lower  left  part  of  the 
maps  with  Ls5r  satisfies  the  assumptions  of  the  Euler-Bernoulli 
theory;  furthermore,  the  assumption  of  essentially  one-dimensional 
phonon  transport  requires  L>lph. 


small  thicknesses,  the  slope  is  fpc.lk^  ()7).k.  As  the  Zener 
model  suggests,  the  position  of  this  peak  can  be  identified 
with  a  characteristic  time  T=flp(!ak,  defining  some  relaxation 
process.  From  the  slope  of  the  peak  position  dependence,  one 
can  further  deduce  scaling  of  r  with  t.  At  large  thicknesses, 
this  characteristic  time  scales  as  t  squared,  which  is  typical 
for  diffusive  processes.  Indeed,  since  case  (a)  presents  a  sce¬ 
nario  of  a  single  Griineisen  constant,  frequencies  of  all 
phonons  are  proportionally  modified  by  strain  and  no  distur¬ 
bances  in  the  local  phonon  distribution  are  created,  although, 
as  a  whole,  the  phonon  distribution  function  corresponds  to  a 
different  temperature,  i.e.,  the  phonon  population  is  either 
heated  or  cooled.  Thus,  the  only  relevant  process  is  the  dy¬ 
namical  spatial  redistribution  of  phonons.  At  large  beam 
thicknesses,  this  spatial  redistribution  is  dominated  by  diffu¬ 
sive  transport  (defining  Fourier  heat  transfer)  and  is  tempo¬ 
rally  slow  with  Tdiff°cf2.  That  is  the  microscopic  phonon  pic¬ 
ture  representation  of  the  classical  thermoelastic  loss 
mechanism  originally  considered  by  Zener.  As  the  beam 
thickness  is  reduced  to  become  comparable  to  the  phonon 
mean  free  path,  the  scaling  of  the  peak  position  changes 
noticeably  to  r«t  to  reflect  a  transition  to  ballistic  phonon 
transfer  across  the  beam  giving  rise  to  the  non-Fourier  heat 
flow,  typical  at  short  times  and  distances. 

Case  (b)  provides  conditions  for  a  clear  observation  of  the 
Akhiezer  effect.  Since  average  y=Piyi+/?2y2=0,  the  local 
phonon  temperature  is  not  changed  and  the  macroscopic  heat 
transfer  does  not  take  place.  On  the  other  hand,  the  local 
phonon  distribution  is  disturbed.  Thus,  the  main  process  is 
the  spatially  local  intergroup  (interbranch)  thermal  relax¬ 
ation,  which  is,  obviously,  defined  by  the  microscopic  char¬ 
acteristic  time  rph.  As  a  result,  the  peak  position  closely  cor¬ 
responds  to  flpeakTph=l  and  shows  no  dependence  on  beam 
thickness  for  a  substantially  thick  resonator.  As  the  sharp 
spatial  inhomogeneities  are  introduced  at  t~  /ph,  the  phonon 
ballistic  transfer  process  becomes,  first,  comparable  in  effi¬ 
ciency  to  interbranch  relaxation  and,  later,  dominates — the 
regime  changes  and  the  peak  position  acquires  dependence 
on  the  mode  frequency. 

Case  (c)  with  y1  =  2yj=0,  72  =  0  demonstrates  the  realiza¬ 
tion  of  a  system  where  both  thermoelastic  and  Akhiezer’s 
mechanisms  are  present  and  of  comparable  strengths.  As  the 
beam  parameters  vary,  their  interplay  results  in  a  complex 
behavior  of  loss  factor,  captured  in  Fig.  1(c).  In  particular, 
two  peaks  are  present  at  large  beam  thicknesses,  which 
merge  and  change  their  character  as  ballistic  transfer  be¬ 
comes  dominant. 

For  a  given  material,  the  pair  of  parameters  7,  and  y2  can 
be  uniquely  specified.  This  can  be  realized  through  careful 
evaluation  of  the  average  values  and  variations  in  the  mode 
specific  Griineisen  constants  available  through  experimen¬ 
tally  tabulated  data  for  different  phonon  branches  or  through 
values  calculated  from  the  second  and  third  order  elastic 
moduli.15,21  It  follows  from  provided  data  that  the  constants 
y1  and  y2  are  typically  highly  unequal  in  magnitude,  with 
say  1 71 1  >  |  y2|,  favoring  the  choice  made  in  case  (c).  It  should 
be  noted  that  our  approach  can  be  easily  generalized  to  in¬ 
corporate  more  than  two  groups  of  phonons,  although  we 
feel  that  this  further  increase  in  complexity  is  unlikely  to 
capture  any  substantially  new  phenomena. 
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To  put  our  numbers  into  perspective,  the  variation  of  Q~x 
over  the  range  of  one  and  a  half  decades  from  its  maximum 
is  explicitly  displayed  in  Fig.  1.  For  the  material  parameters 
chosen,  the  peak  values  of  loss  in  both  the  thermoelastic 
(diffusive)  and  Akhiezer  regimes  correspond  to  Q  —  7.5 
X  104  while  in  the  ballistic  regime,  with  absolute  maximum 
in  the  lower  right  corner,  the  peak  loss  has  a  Q  value  of 
~2.2  X  104. 

Overall,  the  following  quantitative  observations  can  be 
made: 

(i)  Thermoelastic  mechanism  dominates  in  thick  beams 
with  large  aspect  ratio  well  described  by  Euler-Bernoulli 
beams.  In  particular,  with  scaling  down  Si  beam,  the  peak  in 
losses  is  achieved  at  subgigahertz  frequencies,  after  that,  the 
quality  factor  improves. 

(ii)  At  about  f=/ph,  the  character  of  the  phonon  transfer 
changes  to  ballistic,  leading  to  a  different  scaling  of  the  peak 
position  on  the  (II,  f)  map;  this  leads  to  saturation  of  losses 
with  further  reduction  of  the  structure. 

(iii)  Akhiezer  effect  is  of  minor  importance  for  structures 
with  large  aspect  ratio  well  described  by  Euler-Bernoulli 
equation.  It  noticeably  contributes  to  losses  only  at  the  onset 
of  the  regime  of  the  ballistic  phonon  transfer  where  both 
mechanisms  merge. 

(iv)  In  quantitative  terms,  at  1-10  GHz  frequencies,  a  Si 
NEMS  resonator  with  L=5t  has  substantially  higher  losses 
than  a  L  =  lOf  resonator. 


V.  SUMMARY  AND  CONCLUSIONS 

In  summary,  during  gigahertz  frequency  operation,  the 
phonon  subsystem  inside  a  flexural  NEMS  beam  evolves 
into  a  unique  nonequilibrium  dynamics  due  to  a  highly  inho¬ 
mogeneous  strain  pattern.  Resulting  nonequilibrium  heat 
generation  and  redistribution  processes  have  been  treated  be¬ 
yond  the  conventional  heat  diffusion  and  local  temperature 


approximation.  An  advanced  theoretical  analysis  was  formu¬ 
lated  and  conducted  in  the  Boltzmann  framework  to  capture 
the  appropriate  phonon  dynamics  inside  flexural  beams.  In 
the  resulting  study,  the  intrinsic  mechanisms  and  related  ana¬ 
lytic  principles,  including  scaling  rules  influencing  NEMS 
degradation,  have  been  developed  and  expressed  in  a  user- 
friendly  form  for  application. 

The  intriguing  question  of  prohibitive  high -Q  operation 
based  on  fundamental  physical  limitations  and  degradation 
was  also  raised  in  this  study.  Overall,  our  analysis,  using  a 
limited  parameter  set,  shows  no  degradation  in  quality  factor 
with  progressive  miniaturization  from  thermoelastic  and 
Akhiezer  effects;  however,  there  appears  to  be  a  marginal 
softening  in  quality  factor  observed  for  the  Akhiezer  compo¬ 
nent  when  I2rph>  1  as  noted  in  Fig.  1.  Yet,  it  is  commonly 
reported  in  numerical  and  experimental  NEMS  studies  that 
reduction  in  resonator  length  leads  to  substantial  reduction  in 
the  quality  factor  (9. 23-24  These  observations  tentatively  sug¬ 
gest  that  with  down  scaling,  the  microresonator  losses  are 
determined,  to  a  large  extent,  by  surface  related  processes;  in 
fact,  a  larger  fraction  of  the  beam  atoms  reside  at  the  surface, 
so  that  the  surface-to-volume  (SI  V)  component  grows  lin¬ 
early  with  beam  miniaturization.  Therefore,  it  is  likely  that 
an  alternative  mechanism  such  as  enhancement  of  the  third 
and  fourth-order  anharmonic  phonon-phonon  interactions  in 
the  presence  of  the  static  strain  in  the  surface  layer  might 
lead  to  the  direct  dependence  on  the  phonon  characteristic 
time  on  the  surface-to-volume  ratio.  Such  an  anhar¬ 

monic  effect  can  be  incorporated  into  the  general  description 
of  the  phonon  flow  of  this  analysis  through  a  set  of  heuristic 
parameters  describing  scattering  and  thermalization  due  to 
surfaces  or  interfaces.  This  will  be  the  natural  extension  of 
the  present  investigation. 
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